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Abstract — Time-delay  estimation  is  important  in  a  wide  range 
of  applications  in  oceanic  engineering.  In  this  paper,  we  present 
a  novel  time-delay  estimation  algorithm  based  on  maximum  like¬ 
lihood  theory  for  the  case  that  the  measurements  are  corrupted 
by  colored  or  nonuniform  zero-mean  Gaussian  noise.  It  turns  out 
that  the  likelihood  function  associated  with  the  problem  is  highly 
oscillatory,  and  we  propose  a  computationally  efficient  technique 
to  maximize  this  function.  Our  algorithm  first  obtains  an  initial 
estimate  based  on  a  smooth  approximate  cost  function,  and  then 
refines  this  estimate  based  on  the  true  cost  function.  Simulation 
results  show  that  our  estimator  outperforms  a  traditional  phase- 
shift  based  estimator,  and  that  the  estimation  error  approaches  the 
Cramer-Rao  bound  (CRB)  when  the  signal-to-noise  ratio  (SNR) 
increases  without  bound. 

Index  Terms — Supercavitation,  acoustic  sensors,  time-delay  esti¬ 
mation,  maximum  likelihood,  phase-shift  method. 


I.  Introduction 

THIS  paper  deals  with  time-delay  estimation  using  acoustic 
sensors.  The  application,  that  has  served  as  the  main  mo¬ 
tivation  for  this  work  is  the  problem  of  measuring  cavity  thick¬ 
ness  for  underwater  supercavitating  vehicles  (however,  the  re¬ 
sults  of  this  paper  can  also  be  applied  directly  to  many  other 
applications  such  as  proximity  barrier  detection  and  fluid  level 
measurement  [1],  [2]).  In  fact,  the  maximum  speed  of  the  sub¬ 
merged  vehicles,  by  using  traditional  underwater  techniques,  is 
limited  up  to  about  80  m/h.  This  limitation,  however,  is  being 
broken  by  a  physical  phenomenon  called  supercavitation,  which 
makes  it  possible  for  underwater  vehicles  to  travel  at  hundreds 
of  miles  per  hour  [3],  [4],  The  core  of  this  technique  is  to  sur¬ 
round  the  moving  object  with  an  envelope  of  gas  so  that  the 
contact  area  between  the  liquid  and  the  vehicle  surface  is  mini¬ 
mized,  and  as  a  consequence  the  viscous  drag  is  reduced  drasti¬ 
cally.  The  feasibility  of  such  vehicles  is  dependent  mainly  on  the 
ability  to  actively  control  the  shape  and  quality  of  the  gaseous 
cavity  when  the  vehicle  maneuvers  at  high  speeds,  and  hence 
an  efficient  cavity  sensing  technique  is  very  important.  Due  to 
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the  hard  acoustic  boundary  at  the  gas/water  interface,  acoustic- 
based  cavity  thickness  measurement  methods  are  expected  to  be 
superior  to  other  possibilities  such  as  electromagnetic  or  laser- 
based  methods. 

Consider  a  signal  transmitted  from  an  ultrasonic  transducer. 
The  emitted  signal  propagates  through  the  transmission  medium 
(for  instance,  air  or  water)  and  is  reflected  by  the  boundary  (or  a 
suitable  reflector).  The  reflected  signal  is  received  by  a  receiving 
transducer  (that  may  be  the  same  as  the  transmitting  transducer). 
The  time-delay  of  the  signal  propagating  from  the  transmitter  to 
the  receiver  is  measured  and  converted  into  a  value  proportional 
to  the  distance. 

Numerous  techniques  including  the  pulse-echo  and  the 
continuous  wave  (CW)  measurement  techniques  have  been 
developed  for  performing  the  time-delay  measurement  [5]— [7] . 
The  so-called  pulse-echo  method  is  based  on  the  estimation 
of  the  time-delay  between  the  moment  of  transmission  and 
the  time-of-arrival  of  the  reflected  wave.  To  measure  a  small 
distance  and  to  get  a  high  accuracy  of  the  estimate,  the  sound 
pulse  must  be  extremely  short.  However,  owing  to  the  high 
quality  factors  of  the  transducers,  most  existing  ultrasound 
transducers  can  only  generate  the  pulses  with  long  ringing 
tails  [8],  Therefore,  when  using  a  conventional  ultrasound 
transducer,  the  pulse-echo  measurement  technique  is  not  useful 
for  close  proximity  measurements  (that  is,  when  the  distance 
between  the  sensor  and  the  reflector  is  less  than  5  cm). 

Another  class  of  distance  measurement  methods  is  based  on 
continuous  sound  wave  (CW)  generation  techniques  and  over¬ 
comes  the  distance  restriction  of  the  conventional  pulse-echo 
method.  These  methods  require  that  the  transducer  generate  ul¬ 
trasound  continuously  and  receive  the  reflected  signal  simulta¬ 
neously.  The  most  widely  used  method  within  the  class  of  CW 
measurement  techniques  is  called  the  impedance  or  amplitude 
method,  which  can  be  analyzed  similar  to  the  Fabry-Perot  prin¬ 
ciple  known  in  optics  [9],  In  the  space  between  the  sensor  and 
the  reflector,  the  emitted  and  reflected  sound  waves  interfere 
with  each  other  and  a  standing  wave  pattern  is  generated.  The  re¬ 
ceived  signal  amplitude  is  maximal  or  minimal  when  the  waves 
are  in  resonance  or  antiresonance,  respectively.  The  distance 
between  the  sensor  and  the  reflector  is  obtained  by  measuring 
the  frequency  difference  between  two  neighboring  resonance 
peaks  or  notches.  Another  CW-based  method,  referred  to  as 
the  phase-shift  method,  is  much  simpler  and  faster  since  it  only 
requires  two  measurements  of  a  phase-shift  on  two  different 
working  frequencies.  This  method  is  widely  applied  in  optics 
[10] — [12],  A  traditional  estimator  for  the  phase-shift  method 
is  described  in  [6]  and  is  based  on  the  fact  that  given  the  two 
different  working  frequencies,  the  difference  between  the  two 


0364-9059/02$  1 7.00  ©  2002  IEEE 


Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1.  REPORT  DATE 

AUG  2001  2.  REPORT  TYPE 

3.  DATES  COVERED 

00-00-2001  to  00-00-2001 

4.  TITLE  AND  SUBTITLE 

Phase-Shift-Based  Time-Delay  Estimators  for  Proximity  Acoustic  Sensors 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Florida, Department  of  Electrical  and  Computer 

Engineering, Gainesville, FL, 32611 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

18.  NUMBER  19a.  NAME  OF 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE 

unclassified  unclassified  unclassified 

10 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


48 


IEEE  JOURNAL  OF  OCEANIC  ENGINEERING,  VOL.  27,  NO.  1,  JANUARY  2002 


phase-shifts  is  proportional  to  the  time-delay.  However,  it  turns 
out  that  this  estimator  is  not  optimal  when  the  signal  is  known  to 
have  a  real-valued  amplitude.  Furthermore,  this  estimator  does 
not  take  the  noise  into  account  in  a  proper  way. 

In  this  paper,  we  consider  the  phase-shift  method  in  the  case 
when  the  amplitude  of  the  signal  is  nonnegative  unknown,  and 
the  measurements  are  corrupted  by  colored  or  nonuniform 
Gaussian  noise.  A  new  time-delay  estimator  based  on  a  max¬ 
imum-likelihood  approach  is  derived.  However,  this  estimator 
requires  the  solution  of  an  optimization  problem  whose  objec¬ 
tive  function  is  highly  oscillatory.  Inspired  by  a  related  work 
in  [13],  we  propose  an  approach  to  obtain  an  initial  time-delay 
estimate  by  assuming  that  the  amplitude  of  the  reflected  signal 
is  complex-valued  and  then  to  refine  the  initial  estimate  based 
on  the  true  cost  function.  A  computationally  efficient  algorithm, 
referred  to  as  the  direct  matching  method  (DMM),  is  presented 
to  obtain  the  refined  estimate.  As  shown  in  this  paper,  our  al¬ 
gorithm  performs  much  better  than  the  estimator  in  [6]  and  the 
root-mean-squared  errors  (RMSEs)  of  the  parameter  estimates 
of  the  new  method  approach  the  corresponding  Cramer-Rao 
bounds  (CRBs)  when  the  signal-to-noise  ratio  (SNR)  increases 
without  bound.  We  also  show  that  for  the  case  in  which  the 
amplitude  is  no  longer  real-valued,  but  complex-valued  with  a 
small  phase  angle  (i.e.,  in  the  case  of  a  small  model  mismatch), 
the  new  method  can  still  perform  better  than  the  traditional 
method  for  a  low  to  moderate  SNR. 

II.  Problem  Formulation 

Assume  that  the  transmitted  ultrasonic  signal  is  a  sinusoid: 

st(t)  -  siii(27t/if  +  <f> i)  (1) 

where  fi  is  the  transmitted  frequency  and  <f> i  is  the  initial  phase. 
When  the  sound  waveform  is  normally  incident  on  the  interface 
between  two  media,  the  wave  is  partly  reflected  and  partly  trans¬ 
mitted  into  the  second  medium.  The  reflection  coefficient  is  de¬ 
fined  as  R  =  p~  /p+,  where  p+  is  the  pressure  of  the  incident 
wave  and  p~  is  the  pressure  of  the  reflected  wave.  This  can  also 
be  expressed  as  R  —  ( Z2  —  Aj ) / (  a_>  +  z± ),  where  the  z±  and 
Z2  are  the  characteristic  impedances  of  the  transmission  media 
[14].  When  the  impedance  of  the  second  medium  is  much  larger 
than  that  of  the  first  medium,  for  instance,  when  the  aerial  sound 
waves  impinge  onto  a  water  surface,  11  becomes  equal  to  one, 
which  implies  that  the  reflected  signal  is  an  exact  replica  of  the 
incident  signal.  Hence,  in  this  case,  the  signal  received  by  the 
transducer  can  be  written  as 

■sr(t)  —  [1  sin[27t/i(f  -t)  +  4> i  +  7t]  (2) 

where  r  is  the  round-trip  time -delay  between  the  sensor  and  the 
reflecting  boundary  and  (3  is  the  amplitude,  which  is  consid¬ 
ered  to  be  a  nonnegative  unknown,  and  -k  is  induced  by  “half 
wavelength  loss”  due  to  the  acoustic  hard  boundary.  Therefore, 
the  phase-shift  between  the  transmitted  and  received  signal  due 
to  the  time-delay  is  27t/it.  In  order  to  measure  this  phase-shift, 
standard  IQ  processing  is  performed  on  the  received  signal.  This 
means  that  sr(t)  is  demodulated  with  the  in-phase  and  quadra¬ 
ture  components  of  st(t),  sin(27t/if  +  <f> i)  and  cos(27t/if  + 


respectively.  The  demodulated  signals  are  passed  through 
low-pass  filters  and  the  real  and  imaginary  parts  of  f3e-l2'K^lT 
are  obtained  by  measuring  the  output  amplitudes  of  the  two 
low-pass  filters. 

Assume  that  N  subsequent  and  independent  measurements 
are  taken  using  the  frequency  fi .  Then  each  measurement  value 
can  be  modeled  as 

xi(n)  —  f3e->2'K^lT  +  ei(n),  n  —  l,2,...,N  (3) 

where  we  assume  that  e±  (n)  is  a  zero-mean  circular  symmetric 
complex  white  (with  respect  to  n)  Gaussian  random  variable 
with  variance  a2.  Assume  further  that  a  corresponding  set  of 
N  independent  measurements  are  obtained  by  transmitting  a 
signal  with  the  known  frequency  /2  (without  loss  of  generality, 
we  assume  /;  <  /2): 

X2{n)  =  f3e3‘2'Kf'2T  +  e2{n),  n  =  1,2,  ...,7V  (4) 

where  e2 (n)  is  Gaussian  noise  with  variance  a\.  We  will  in  gen¬ 
eral  assume  that  the  noise  is  nonuniform  so  that  a2  ^  a2.  This 
modeling  accounts  for  the  possibility  that  the  noise  level  is  dif¬ 
ferent  at  different  frequencies,  a  situation  that  can,  for  example, 
arise  from  hardware  artifacts.  Furthermore,  we  will  allow  for 
the  possibility  that  ei(n)  and  e2  in)  are  correlated,  a  situation 
that  could  arise,  for  instance,  if  the  two  measurements  (at  time 
n)  on  frequency  fi  and  /2  are  taken  simultaneously  using  two 
transducers.  We  will  see  that  this  data  model  leads  not  only  to 
a  robust  estimation  algorithm,  but  is  also  mathematically  effi¬ 
cient. 

The  data  model  can  be  expressed  using  matrix  notation  as 
follows.  Define  the  vector  of  measurements 

x(n)  =  [£i(n)  X2 (n)]T  6  C2xl  (5) 

where  ( •  )T  denotes  the  transpose,  let 

e(n)  =  [ei(n)  e2(n)]T  6  C2xl  (6) 

be  a  vector  containing  the  noise  samples  and  let 

b(r)  =  [eJ,27r^lT  ej2vf2T]T.  (7) 

Then  the  data  model  in  (3)  and  (4)  can  be  written  in  a  matrix 
form  as 


x(n)  =  fib  +  e(n),  n  —  1, 2, . . . ,  N  (8) 

where  e(n)^'_1  are  independently  and  identically  distributed 
zero-mean  circularly  symmetric  complex  Gaussian  random  vec¬ 
tors  with  an  unknown  positive  definite  covariance  matrix 


Q  = 


of  poqoA 
A*  CT1CT2  O-f 


(9) 


and  here  the  complex  correlation  coefficient  p  is  defined  as 


E{e\e2) 

P  = - 

oqoA 


(10) 


Hereafter  E(x )  is  the  expected  value  of  x  and  ( •  )*  denotes  the 
complex  conjugate.  The  problem  of  interest  in  our  paper  is  to 
estimate  {r,  (3}  from  x(n),  n  —  1,2, ...  ,N. 
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III.  Time  Delay  Estimation 

A.  Maximum  Likelihood  Estimator 

The  log-likelihood  function  of  x(n)  is  proportional  to  (within 
an  additive  constant): 

f  i  N 

C  =  -  ln[det(Q)]  -  tr  <  Q  1  —  ^[x(n) 

V  n= 1 

-  /?b][x(n)  -  ib]7i  |  (11) 

where  det(  • )  denotes  the  determinant  of  a  matrix,  [  •  \H  denotes 
the  conjugate  transpose,  and  tr{-}  stands  for  the  trace  of  the 
matrix. 

It  is  easy  to  show  that  the  matrix  Q  that  maximizes  (11)  is 
(see,  e.g.,  [15],  [16]) 

1  N 

Q  =  —  ^[x(n)  -  /3b][x(n)  -  (3b\H .  (12) 


Before  we  proceed,  let  us  remark  on  the  following  two  facts. 
First,  the  data  model  in  (8),  the  related  optimization  problem 
of  maximizing  (11)  and  its  solution  are  similar  to  those  ob¬ 
tained  for  the  amplitude  and  phase  estimation  ( APES )  filter  [  1 8  ] . 
Second,  it  is  evident  from  (20)  that  the  exact  ML  estimator  of 
(t.  (3)  reduces  to  a  simple  nonlinear  least  square  (NLS)  cost 
function  of  pre-whit ened  data.  Note,  however,  that  Q  in  (19) 
is  not  the  ML  estimate  of  Q. 

We  will  consider  three  different  cases:  [3  is  complex-valued,  [3 
is  real-valued  and  (3  is  nonnegative,  respectively.  Consider  first 
the  case  where  (3  is  an  unknown  complex-valued  scalar.  Then 
minimizing  C2(r,  (3)  with  respect  to  t  and  (3  yields  the  estimates 
of  f  and  / 3 ,  respectively,  as 


A  /  \  A  bffQ  xx  - 

rCl  =  argmaxcijr )  =  argmax  A 

T  T  bffQ“1b 

a  _  bgQ-1x 

/ci  bffQ_1b  ' 


(21) 

(22) 


Insertion  of  (12)  into  (11)  shows  that  maximizing  (1 1)  is  equiv¬ 
alent  to  minimizing 

Ci  =  det  ^  ^j[x(n)  _  ^b][x(n)  -  /?b]ff^  =  det(G) 

(13) 

where 

1  N 

G  =  —  -  /3b][x(n)  -  (3hf .  (14) 

n= 1 

Let 

1  N 

x=]yXXn)’  (15) 

n=  1 

and 

1  ^ 

rm  =  ^^2  x(n)x(n)ff  •  (l6) 

n=  1 

Then  G  can  be  written  as 

G  =  Rxx  -  [3bxH  -  f3*xbH  +  \f3\2bbH 

=  [/  lb  -  x]  [/lb  -  x]7i  +  R.„;  -  xxJi.  (17) 

With  this  notation,  the  cost  function  C\  in  (13)  can  be  written  as 

Ci  =  det[RM  -  xxff  +  (Jib  -  x)(/?b  -  x)ff] 

=  det(Q)  det  [I  +  Q_1(/'ib  —  x)(/?b  —  x)ff] 

=  det(Q)[l  +  (Jib  -  x)ffQ-1(/?b  -  x)]  (18) 

where 

Q  =  RM-xxff  (19) 


and  where  we  have  used  the  fact  that  det  (I  +  AB)  =  det  (I  + 
BA)  whenever  the  dimensions  of  A  and  B  are  conformable 
[17].  Hence,  maximizing  C  is  equivalent  to  minimizing 

C2(t,  (3)  -  \J3b{r )  -  x]iAQ_1[/'ib(r )  -  x]  (20) 

which  is  a  highly  nonlinear  optimization  problem. 


Next,  consider  the  case  of  a  real-valued  (3.  Minimizing  C2(t ,  (3) 
with  respect  to  t  and  (3  yields 

,  a  Re2(bffQ_1x) 

tCo  —  argmax  cAt)  —  argmax - * - 

r  r  bffQ-1b 

3  _  Re(bffQ_1x) 

°2  ~  bffQ_1b 

T-TC2 

where  Re(at)  is  the  real  part  of  x.  Finally,  consider  the  case  when 
[3  is  a  nonnegative  unknown.  Minimizing  C2(r ,  (3)  for  a  fixed  r 
with  respect  to  [3  yields 


(23) 

(24) 


_  ReQ^Q^x)  / ReQj^Q^xA 

b^Q-ib  U  \  b^Q-ib  J 


(25) 


where  u(x)  is  the  unit  step  function: 


u(x)  — 


1,  if  x  >  0 
0,  if  x  <  0. 


Insertion  of  (25)  into  (20)  gives 


(26) 


fC3  =  argmax  03(7-) 

r 


A 

=  argmax 


Re2(bffQ_1x) 

- ^ - u 

b^Q-ib 


f  ReO^Q^xA 
\  bffQ_1b  )  ' 


(27) 


Once  we  have  obtained  tC3  ,  the  corresponding  Hc,  is  obtained 
from  (25)  as  f%3  —  /5(fC3). 

According  to  the  parsimony  principle  [19],  if  (3  is  real-valued, 
the  estimates  obtained  by  maximizing  c2(r)  and  c^(r)  should 
be  more  accurate  than  those  obtained  by  maximizing  ci  (t)  due 
to  the  different  constraints  on  [3.  Although  (23)  and  (27)  are 
simple  one-dimensional  search  problems,  the  objective  func¬ 
tions  c2(t)  and  cs(t)  are  highly  oscillatory  and  have  numerous 
closely  spaced  local  maxima  which  make  it  very  difficult  to  find 
the  global  maximum.  In  the  section  below,  we  present  a  fast  al¬ 
gorithm  to  cope  with  these  optimization  problems. 


50 


IEEE  JOURNAL  OF  OCEANIC  ENGINEERING,  VOL.  27,  NO.  1,  JANUARY  2002 


B.  Estimation  Algorithm 

To  solve  the  optimization  problems  in  previous  section,  we 
first  introduce  the  following  notation.  Let 


a  r  71  v 

(28) 

L v  72  J 

-  ai'g('y) 

(29) 

1  N 

*  =  1,2 

(30) 

n— 1 

,\ 

=  arg(Ti),  i  =  1,2 

(31) 

where  Q  is  defined  in  (19)  and  aig(x)  is  the  phase  of  x. 

1)  Maximization  of  Ci(r):  Consider  first  the  optimization 
problem  for  the  cost  function  Ci(r),  which  after  some  straight¬ 
forward  algebra  can  be  written  as 

,  ,  A+\B\cos(iI’  +  i/’b) 

1  7i  +  72  +  2H  cos (V’  +  Vv) 

where 

A^|a7|2(7f  +  H2)+I^|2(7l  +  H2) 

+  2Re[TiT2(7i  +  72)^*]  (33) 

B  =  2[;y(|Ti|27i  +  |72|272)  +  77*271 72  +  T^2^2](34) 

4’b  =  ai'g  (B),  (35) 

V’  -  27 r(/2  -  fi)T.  (36) 

Differentiating  <7  (r)  with  respect  to  7  and  setting  the  derivative 
to  zero  yields 

\B\  sin(V’  +  V’b)[7i  +  72  +  2\u\  cos(f  +  Vv)] 

-[A  +  \B\ cos(V’  +  4’bWW\  sin(V’  +  Vv)]  =  0.  (37) 


p  -  | i? | (71  +  72)  cost’s)  -  2|;/|Acos(Vv)  (38) 

q  =  |S|(7i  +  72)  sin(V’s)  -  2|;y|Asin(Vv)  (39) 

4/  =  arg(p  +  jq).  (40) 

Then  the  solution  of  (37)  is 

-  .  (  2|;y||i?|  siii(V’s  -  Vv)  \  T  ,,,, 

*  = arcsm  { - —  j  “  *■  (41) 

If  it  is  known  a  priori  that  rmax  <  1  /(/2  —  /i),  the  estimate  of 
r  is 

if^° 

T<*  =  \  V-  +  2  7T  .  .  (42) 

Note  that  in  the  special  case  that  it  is  known  that  Q  =  o2 1  where 
I  is  the  identity  matrix  and  o2  =  o2  =  a2,  the  estimator  in  (42) 
can  be  simplified  to 

f  V’a2  ~  V’ai  if  7-  -7-  >0 

A  -  J  2n{h-h)  ’  11  *2  (.~ 

01  1  (i-’xi  —  iA,)  +  2tt  ,  .  „ 

l  "277/2-/1)  ’  lfV’®2  V’*l  <  0 

which  is  the  same  expression  as  for  the  traditional  phase-shift 
method  in  [6], 


2)  DMM:  Although  a  closed-form  solution  to  maximizing 
ci(t)  with  respect  to  r  can  be  found,  it  is  less  accurate  than 
the  solution  obtained  by  maximizing  c2(r)  and  c${r),  if  fi  is 
real-valued.  However,  it  is  difficult  to  maximize  c2(r)  and  03(7") 
directly  due  to  the  fact  that  both  of  them  are  highly  oscillatory 
with  Cz(t)  being  even  nonsmooth. 

In  this  paper,  we  propose  to  maximize  the  cost  function  <7  (r ) 
by  first  obtaining  an  initial  estimate  fCl ,  and  then  refine  it  based 
on  c2(t)  or  c^(t)  respectively.  Our  algorithm,  which  doesn’t 
require  a  search  through  the  parameter  space  of  r,  is  described 
below  and  will  be  called  the  DMM  algorithm. 

Rewrite  the  cost  function  c2(r)  as 

c2{j )  =  <7i(t)  •  <72(t )  (44) 

where 

9i(j)  =  (bffQ_1b)_1 

=  [71  +  72  +  2|z/|  cos(27 r(/2  -  /i)r  +  Vv)]_1  (45) 

and 

72(7- )  =  Re(bffQ_1x) 

=  7i  |*i  |  cos(27t/ir  -  ) 

+  72|37|  cos(27r/2r  -  V’s2) 

+  |z/|  |7i  |  cos(2t r/2r  -  V’®!  +  Vv) 

+  |//||72|  cos(2tt/it  -  V’s2  -  Vv)-  (46) 

It  can  be  noted  that  the  denominator  of  gi  (r)  contains  a  sinusoid 
with  (low)  frequency  /2  —  fi  and  that  g2(j)  is  highly  oscilla¬ 
tory  with  an  oscillation  frequency  somewhere  between  fi  and 
/2.  Therefore,  the  cost  function  c2(r)  is  highly  oscillatory.  Dif¬ 
ferentiating  c2(t)  with  respect  to  r  yields 

4(r)  =  9i(r)  ■  9l{r)  +  2gi{r)  ■  g2(r)  ■  gf2{r).  (47) 

Next,  we  expand  c'At)  around  TCl  to  a  first-order  approxima¬ 
tion  and  find  the  solution  to  the  equation  c'A  t)  —  0,  which  cor¬ 
responds  to  the  local  maximum  of  c2(r)  around  fCl.  A  Taylor 
series  expansion  of  <7i(r)  around  fCl  yields  to  a  second-order 
approximation 

9i (t)  wa0  +  ai(r-  fCl)  +  a2(r  -  fCl)2  (48) 


where 


«0  —  9i(.t)\t—tc1 

-  [71  +  72  +  2|//|  cos(2t r(/2  -  fi)fCl  +  Vv)]_1 


a  dgiij) 

al  ~~  - o - 

or 


'■  4tt(/2  -  /i)H  sin(27r(/2  -  /i)fCl  +  Vv)[7i 
+  72  +  2|i/|  cos(2t r(/2  -  /i)fCl  +  Vv)]-2 


(49) 


(50) 


A  192^i(r) 

fl2=  2^^r  f 

T— Tc^ 

=  4tt2(/2  -  /i)2H{cos(27r(/2  -  /i)fCl  +  Vv)[7i 
+  72  +  2|z/|  cos(2t r(/2  -  /i)TCl  +  Vv)]-2 
+  4|z/|  sin2(27r(/2  -  /i)fCl  +  Vv)[7i  +  72 
+  2|z/|  cos(2t r(/2  -  /i)fCl  +  Vv)]-3}-  (51) 
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Since  fi  is  nonnegative,  we  have  the  following  approximation 
for  high  SNR: 


2vr/ifCl  ~  2kjTY  +  4% 

*  =  1,2 

(52) 

where 

kj.  — 

2vr/ifCl  -  V’*i 

\i  *  =  1,2 

(53) 

2tt 

with  [  •  J  denoting  rounding  to  the  nearest  integer.  Note  that  the 
information  that  (3  >  0  is  used  explicitly  here.  In  a  small  range 
around  fCl,  (52)  yields  the  following  approximation: 

cos[2t r/i(T  -  fCl)]  ps  cos(2t t/,:t  - 

:%xi  -  PS  1, 

*  =  1,2 

(54) 

and 

sin[27T fi(r  -  TCl)\  ps  sin(27r/,T  - 

4’xi  -  2Atf7T ) 

PS  27T fiT  -  4’Si 

—  2/777  i  —  1,  2. 

(55) 

Then  uf  j)  and  gif.T)  can  be  expanded  in  their  first-order  ap¬ 
proximations  around  fCl  as 

92{t)  PS  60  +  biT 

(56) 

and 

gift)  ps  do  +  d\T 

(57) 

where 


bo  —  |*i  |  [71  +  M  cos(Y’)  -  \v\  sin(V’)(V’*2  +  2fc27t)] 

+  1^2!  [72  +  M  COS  (V~)  +  \v\  sill  (VOW'S !  +  2fci7r)] 

(58) 

61  =  2-k\v\  sin(V’)(/2|^i|  -  /i|*2|)  (59) 

do  -  27r|;y|sm(V’)(/2^i|  -  /ll^l) 

+  27r/i(V’*1  +  2fci7r)[|;yp2|  cos(V’)  +  7i4’i|] 

+  2tt /2(V’s2  +  2fc27r)[|;ypi|  cos(V’)  +  72^2!]  (60) 

4  =  -47t2|fi|  [|u|/|  cos(Y’)  +  71/f] 

-  47t2|7:2|  [|^|/f  cos(V’)  +  72/2]  (61) 

where 

V"  =  V’®!  -  4’x2  -  4'w  (62) 

Therefore,  around  fCl ,  c2(t)  can  be  approximately  expressed  as 


and 


doi  —  do  +  d\TCl . 


(66) 


As  a  final  remark,  we  note  that  if  it  is  known  a  priori  that  Q  = 
a2 1,  i.e.,  71  =  72  and  1 /  =  0,  and  if  the  SNR  is  high,  i.e., 
N  .r2|,  (64)  can  be  simplified  to 


tdmm|q=<t2i 


fl  (  2fcl7T  +  j’Sl  \ 
ft  +  /i  v  27T/!  ; 

,  ft  (  2fc27T  +  V’x2  \ 

ft  +  ft\  2t/2  j- 


(67) 


Equation  (67)  can  be  interpreted  as  a  weighted  summation  of 
the  two  estimates  of  r  obtained  via  fi  and  f  >  separately. 

The  steps  of  the  DMM  algorithm  are  summarized  as  follows: 

1)  Estimate  Q  via  (19).  Invert  Q  to  obtain  71, 72  and  v  [see 
(28)]. 

2)  Obtain  an  initial  estimate  of  fCl  via  (42). 

3)  Estimate  ki  and  k 2  via  (53). 

4)  Compute  a®,  ai,  a2,  bo,  bi,  do,  d±  via  (49)-(51)  and 
(58)-(62). 

5)  Obtain  the  final  estimate  tdmm  via  (64). 

Once  the  final  estimate  f  is  obtained,  computing  the  estimate 
fi  is  straightforward  via  (24). 

3 )  Analysis  of  Model  Errors:  In  practical  measurement  sce¬ 
narios,  due  to  the  complex  measurement  environment  and  in 
particular  the  nonideal  rigid  wall  reflection  and  inhomogeneities 
of  the  transmission  media,  it  may  not  be  entirely  correct  to 
model  fi  as  a  nonnegative  number.  In  fact,  it  can  be  argued 
that  fi  should  be  modeled  as  a  complex-valued  number  fi  — 
\ii\cJ,'  :  where  iff s  is  a  small  phase  angle.  In  this  case,  the  esti¬ 
mate  tdmm-  which  is  obtained  by  assuming  (}  is  nonnegative, 
will  be  biased  due  to  the  model  mismatch.  It  can  be  shown 
that  for  the  special  case  Q  =  a2 1  in  (67),  the  bias  and  the 
mean-squared  errors  (MSE)  of  the  estimate  tdmm  are 


Bias(foMM)  =  £(?dmm  ~  TIT)  — 

and 


/1  +  /2 

Mft  +  ft) 


4f)  (68) 


MSE(tdmm)  =  £[(tdmm  -  r)2|r]  = 


m\[i\v-{ft  +  ft) 


/1  +  /2 
2vr(/i2  +  fi) 


1  2 


if) 


(69) 


4 (t)  ~  [ai  +  2 a2(r  -  fCl)](60  +  6it)2 

+2[a0  +  «i(t  -  fCl)](60  +  61  t )(d0  +  d\T ).  (63) 

Simplifying  (63),  neglecting  the  second-order  term  of  t  —  fCl , 
and  setting  it  to  zero,  we  obtain  the  refined  solution,  as  shown 
in  (64)  at  the  bottom  of  the  page,  where 

boi  -b0  +  biTCl  (65) 


It  is  clear  from  the  above  equations  that  both  Bias(fDMM)  and 
MSE(tdmm)  increase  as  if 3  increases. 

IV.  Numerical  Examples 

In  this  section,  we  present  several  computer  simulations  to  il¬ 
lustrate  the  performance  of  our  proposed  algorithm.  The  perfor¬ 
mance  of  our  algorithm  is  compared  to  the  Cramer-Rao  bound 


T DMM  —  TCi 


aib201  +  2aoboidoi 

2  ((12^01  4”  aiMoi  +  aobidoi  +  aoboidi  +  oi&oi^oi) 


(64) 
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1  1.5  2  2.5 

Time  delay  x  (Second) 


X  10 


(a) 


(b) 


(c) 


Fig.  1.  Comparison  of  the  cost  functions  f  i(  r),  c2(r)  and  c3(r)  for  different  assumptions  on  the  signal  amplitude  i.  No  noise  is  present.  True  time-delay  is 
1.67  X  10-4  s.  (a),  (b)  f)  is  nonnegative,  three  cost  functions  share  same  global  maximum,  (c)  >  is  complex-valued,  the  global  maximum  of  c3(r)  deviates  away 
from  the  true  location  of  r . 


(CRB),  which  gives  the  minimum  attainable  variance  for  any 
unbiased  estimator.  The  CRB  for  our  problem  is  derived  in  Ap¬ 
pendix  A.  In  the  simulations  below,  we  use  N  =  10,  [3  = 
0.5, /i  =  60  kHz,  f2  —  63  kHz,  r  =  1.67  X  10-4  s  (corre¬ 
sponding  to  0.028  m  if  the  sound  speed  is  340  m/s).  The  RMSE 
of  each  estimate  is  obtained  by  500  Monte  Carlo  trials  and  the 
SNR  is  defined  as  2 [32  /(a2  +  aI)- 

Fig.  1  illustrates  the  cost  functions  Ci(r),  c2(t)  and  C3 (t).  As 
shown  in  Fig.  1(a)  and  (b),  Ci(t),  c2{t)  and  03(7")  share  the  same 
global  maximum  if  fi  is  nonnegative.  Note  that  both  c2(t)  and 
c2(t )  are  highly  oscillatory,  and  also  that  03(7")  has  many  non- 
differentiable  points.  On  the  other  hand,  Ci(r)  is  the  envelope 
of  c2(t )  and  C3 (r)  and  is  in  fact  quite  smooth.  It  is  expected  that 
maximizing  c2(t)  or  03(7")  can  yield  a  much  more  accurate  es¬ 
timate  than  maximizing  ci  (r),  due  to  the  sharper  peaks  of  c2(t) 
or  c2{t).  Moreover,  since  c2{t)  has  more  closely  spaced  local 
maxima  than  c2(t),  it  is  intuitively  expected  that  c2{t)  is  more 
sensitive  to  noise  in  a  certain  SNR  range  than  C3 (r )  is.  If  (3  is 
complex-valued,  as  shown  in  Fig.  1(c),  the  global  maxima  of 
c2 (t )  and  03(7")  deviate  from  the  true  location  of  r  (for  sim¬ 


plicity,  only  c2(t)  is  shown  in  this  figure).  Hence,  the  estimates 
based  on  c2(t)  or  03(7")  will  be  biased  due  to  the  model  mis¬ 
match. 

Fig.  2(a)  illustrates  the  RMSEs  of  the  different  methods  for 
the  estimates  of  r.  The  noise  covariance  matrix  used  here  is  Q  = 
a2 1  and  five  different  methods  are  compared.  The  first  is  called 
“Cl”  and  uses  the  closed-form  solution  in  (42)  to  maximize  cost 
function  Ci(r).  The  second  and  the  third  methods  are  based  on 
maximizing  c2(t)  and  03(7")  by  a  1-D  search  method,  and  are 
referred  to  as  “C2”  and  “C3,”  respectively.  The  1-D  search  is 
performed  in  two  steps,  with  an  initial  estimate  via  (42),  fol¬ 
lowed  by  a  fine  search  using  the  / mill  function  in  MATLAB. 
The  fourth  is  “C3  (White),”  which  is  a  special  case  of  C3,  where 
the  Q  is  set  to  an  identity  matrix,  that  is,  we  assume  the  noise 
is  white.  The  last  method  is  the  DMM  given  by  (64)  [(67)  is  not 
used  in  any  of  our  simulations].  It  can  be  noted  that  Cl  cannot 
approach  the  real  CRB  even  when  the  SNR  goes  to  infinity. 
However,  both  C2,  C3,  and  DMM  approach  the  real  CRB  when 
the  SNR  increases.  Note  that  the  threshold  effect  for  C3  and 
DMM  occurs  earlier  than  for  C2,  which  is  intuitively  expected 
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(a) 


(b) 


Fig.  2.  Comparison  of  the  RMSEs  and  CRBs  for  estimating  (a)  r ;  (b)  (3  when  the  signal  amplitude  is  assumed  to  be  complex-valued  (Cl),  real-valued  (C2),  and 
nonnegative  (C3).  The  noise  covariance  matrix  is  Q  =  <r2I. 


Fig.  3.  Comparison  of  the  RMSEs  and  CRBs  for  estimating  r  when  the 
signal  amplitude  is  assumed  to  be  complex-valued  (Cl),  real-valued  (C2),  and 
nonnegative  (C3).  The  noises  of  the  two  channels  are  correlated  with  each 
other  and  have  different  variances. 

since  both  C3  and  DMM  use  the  information  that  [}  is  nonneg¬ 
ative.  Note  also  that  the  performances  of  C3,  C3(White)  and 
DMM  are  almost  the  same,  despite  the  fact  that  C3  and  DMM 
do  not  take  advantage  of  the  prior  knowledge  of  the  noise.  In 
other  words,  the  robustness  against  the  colored  or  nonuniform 
noise  offered  by  our  data  model  and  algorithm  comes  at  a  neg¬ 
ligible  cost  (or  rather  a  negligible  performance  loss  in  the  case 
Q  =  a2 1).  The  reason  for  this  is  related  to  the  fact  that  our 
estimation  problem  is  decoupled  in  the  sense  that  the  Fisher  in¬ 
formation  matrix  is  block  diagonal  (cf.  Appendix  A).  Fig.  2(b) 
shows  the  corresponding  estimates  of  fj. 

Fig.  3  illustrates  the  RMSEs  of  the  different  methods  for 
the  estimates  of  r  when  the  noise  covariance  matrix  is  Q  = 
[—0.93 Iwj  _0-9^1{72j]  where  af/aj  =  2.  Note  that  the  per¬ 
formance  of  DMMis  close  to  that  of  C3  and  approaches  the  real 
CRB  as  the  SNR  increases.  If  the  information  of  Q  is  not  used. 


as  shown  by  C3  (white),  there  is  a  gap  of  8  dB  from  the  real 
CRB,  no  matter  how  high  the  SNR  is. 

As  a  further  illustration  of  the  robustness  of  our  algorithm. 
Fig.  4(a)  shows  the  performance  of  the  C3  and  DMM  as  a  func¬ 
tion  of  the  absolute  value  of  the  correlation  coefficient  p  by 
fixing  the  SNR  =  35  dB,  arg(p)  =  7t/2  and  a2  —  a2.  The 
figure  shows  that  both  C3  and  DMM  are  close  to  the  real  CRB, 
however,  C3  (White)  deviates  from  the  real  CRB  when  \p\  ap¬ 
proaches  unity  (i.e.,  when  the  noises  in  the  two  channels  become 
highly  correlated).  Fig.  4(b)  shows  the  performance  of  the  r  es¬ 
timators  as  a  function  of  the  ratio  of  the  two  variances  o2/o2 
by  fixing  SNR  =  35  dB  and  p  —  0.  It  can  be  noted  that  C3 
and  DMM  can  still  achieve  the  real  CRB,  whereas  C3  (white) 
deviates  from  the  real  CRB  when  the  ratio  increases.  The  inter¬ 
pretation  of  this  is  that  our  estimator  has  ability  to  estimate  the 
noise  level  on  the  two  different  frequencies,  and  to  combine  the 
measurements  optimally  taking  the  different  noise  levels  into 
account. 

From  Figs.  3  and  4,  we  can  see  that  our  new  estimator  has  the 
ability  to  adapt  to  an  unknown  noise  model.  Both  C2,  C3,  and 
DMM  can  achieve  the  real  CRB  as  long  as  the  SNR  is  above 
a  certain  threshold.  However,  DMM  is  computationally  more 
efficient  than  C2  and  C3  since  it  does  not  require  any  search. 

Finally,  in  Fig.  5  we  illustrate  the  performance  of  DMM  when 
rl  is  complex-valued  with  a  small  phase  angle.  The  noise  covari¬ 
ance  matrix  is  Q  =  o2 1.  Due  to  the  model  mismatch,  our  esti¬ 
mates  are  biased,  and  can  therefore  not  approach  the  real  CRBs 
as  the  SNR  increases.  However,  if  the  phase  angle  i/’/t  is  small 
and  the  SNR  is  low  to  moderate,  which  means  that  the  estima¬ 
tion  error  caused  by  the  noise  is  larger  than  the  constant  bias 
due  to  the  model  mismatch,  DMM  still  performs  better  than  the 
traditional  method.  When  the  SNR  increases,  the  constant  bias 
caused  by  the  model  mismatch  becomes  dominant  compared  to 
the  error  caused  by  the  noise,  and  the  traditional  method  out¬ 
performs  DMM.  Note  that  the  RMSE  of  fi  for  the  traditional 
method  is  lower  than  the  corresponding  CRBs  at  the  low  SNRs, 
which  is  due  to  the  fact  that  it  is  a  biased  estimator. 
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(a)  (b) 


Fig.  4.  Comparison  of  the  RMSEs  and  CRBs  for  estimating  r .  (a)  As  a  function  of  the  absolute  value  of  p  with  arg(p)  =  7r/2 .  (b)  As  a  function  of  <r\ /<j\  with 
p  =  0.  The  signal  amplitude  is  assumed  to  be  nonnegative  and  SNR  =  35  dB. 


SNR  (in  dB)  SNR  (in  dB) 


(a) 


(b) 


Fig.  5.  Comparison  of  the  RMSEs  and  CRBs  by  using  DMM  to  estimate  (a)  r  and  (b)  fi.  Here  the  phase  angle  of  fi  is  0  ' .  I11,  2  ,  4  ,  respectively.  The  noise 
covariance  matrix  is  Q  =  <r2I. 


V.  Conclusion 

In  this  paper,  we  have  presented  a  new  time -delay  estimator 
for  the  phase-shift  method.  We  have  considered  the  case  when 
the  amplitude  of  the  signal  is  nonnegative  unknown  and  the  mea¬ 
surements  are  corrupted  by  Gaussian  noise  with  unknown  co- 
variance.  To  deal  with  the  induced  optimization  problem,  whose 
cost  function  is  highly  oscillatory,  we  propose  a  method  called 
DMM  that  first  finds  an  initial  estimate  of  the  unknown  param¬ 
eter  by  maximizing  a  smoother  cost  function,  and  then  refines 
this  estimate  based  on  the  true  cost  function.  The  RMSEs  of  the 
estimates  obtained  by  using  our  new  algorithm  approach  the  cor¬ 
responding  CRBs  as  the  SNR  increases.  Since  DMM  does  not 
require  any  search  over  the  parameter  space,  it  is  suitable  for 
a  practical  implementation.  We  also  show  that  in  the  case  of  a 
model  mismatch,  if  the  phase  angle  of  the  reflection  coefficient 
is  small,  the  DMM  can  still  perform  better  than  the  traditional 
method  for  low  to  moderate  SNR.  Although  the  motivation  for 
our  work  has  been  cavity  thickness  measurements  for  supercav¬ 


itation  systems,  the  results  of  this  paper  are  applicable  to  many 
other  proximity  distance  measurement  applications  as  well. 

Appendix  A 
Derivation  of  CRBs 

We  outline  below  the  derivation  of  CRBs  for  the  parameter 
estimates  of  the  following  data  model: 

x(n)  =  fib  +  e(n),  n  —  1, 2, . . . ,  N  (70) 
where  b  =  [e-,2ir^lT  e427r/2T]X. 

In  (70)  the  additive  noise  vectors  e(n),  n  —  1,  2, . . . ,  N  are 
assumed  to  be  independent  zero-mean  Gaussian  random  vec¬ 
tors  with  an  unknown  covariance  matrix  Q.  The  data  samples 
x(n),  n  —  1,2,. . ,  ,1V  are  independent  and  the  unknown  pa¬ 
rameters  in  the  estimation  problem  are  the  real  and  imaginary 
parts  of  fi  (for  complex-valued  fi)  or  simply  the  amplitude  of  fi 
(for  real-valued  /  (),  the  time-delay  r  and  the  unknown  elements 
of  Q.  Let  7/  be  a  vector  containing  those  parameters.  The  CRB 
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inequality  states  that  provided  the  estimates  are  unbiased,  the 
covariance  matrix  of  the  estimated  parameter  vectors  //  is  lower 
bounded  by  the  inverse  of  the  Fisher  information  matrix  FIM-1 . 
The  FIM  can  be  computed  by  the  extended  Slepian-Bang’s  for¬ 
mula  [17]: 

{FIM},,-  =  (Vtr(Q_1Q-Q_1Q}) 

+2N  Re[(/ib) Q_1  (/3b)}]  (71) 

(here  X]  denotes  the  derivative  of  X  with  respect  to  the  zth  un¬ 
known  parameter).  Note  that  the  FIM  is  a  block  diagonal  matrix 
since  Q  does  not  depend  on  the  parameters  in  pb(r),  and  /9b  (r) 
does  not  depend  on  the  elements  in  Q.  Hence,  the  CRB  of  the 
estimates  of  the  delay  and  the  amplitude  can  be  determined  from 
the  second  term  of  the  right  side  of  (71). 

We  first  derive  the  CRB  s  for  the  case  of  a  complex-valued  fj. 
Let 


3S3 

o 

II 

¥ 

1 — 1 

(72) 

where  IiiiLc)  denotes  the  imaginary  part  of  x,  be  a  vector  of  the 
signal  parameters  (i.e.,  the  parameters  in  the  upper  left  block  of 

the  FIM).  Define 

.  /  A  dr-  u 

Ml  0Re(/3) 

(73) 

/  a  dp,  ,u 

A'-  ~~  dlm(p)  ~  J 

(74) 

/  A  dp,  _  \fieP^hr~ 

F'3  qt  J9ej2nf2T 

(75) 

and  let 


Let 


F  =  K  p'2\.  (83) 

Then 


CRB(r/sr)  =  [27VRe(FffQ-1F)]-1.  (84) 


If  a2  =  a2  =  a2  and  p  —  0,  simplifying  (84)  yields 


GRB"  =  4JV 


and 


CRBr  = 


mpv-{f2  +  f2y 


(85) 

(86) 


In  general,  the  frequency  difference  is  determined  by  the 
maximum  distance  we  want  to  measure.  Therefore,  once 
l/i  —  /2I  is  fixed,  we  can  choose  large  fi  and  f2  to  decrease 
the  real  CRB  of  r.  Note  however  that  the  lower  real  CRB  we 
design,  the  higher  SNR  is  required  to  approach  it. 

Note  that  in  the  case  when  Q  =  a2 1,  according  to  (79)  and 
(86),  the  ratio  of  the  CRBs  of  r  for  complex-valued  [}  and  for 
real-valued  [i  is 


CRBcomplex  =  2  (f2  +  /|) 
CRBreal  (h  -  f2f  ’ 


This  ratio  is  directly  related  to  the  ratio  of  the  oscillation  fre¬ 
quencies  of  ci(t)  and  c2(r),  respectively,  and  therefore  to  the 
sharpness  of  the  peaks  of  Ci(r)  and  c2(r).  This  confirms  the 
heuristic  reasoning  in  Section  IV  which  concludes  that  the  cost 
function  with  the  sharpest  peak  should  give  the  most  accurate 
estimate. 


F  —  [f-1  f4  f4]-  (76) 


Then 


CRB(r/sc)  =  [2IVRe(FffQ_1F)]_1  (77) 

which  is  easily  evaluated  numerically.  If  a2  =  a2  =  a2  and 
p  —  0,  simplifying  (77)  yields 


CRBd  —  CRBRe(d)  +  CRBim(d) 

=  £(1  +  2  JtUL 

4  (fi  -  f2)2 


and 


CRBT  = 


4N\f3\2ir2(fi  -  f2 )2 


(78) 

(79) 


Clearly,  for  maximal  accuracy,  f2  —  fi  should  be  as  large  as 
possible.  Note  however  that  f2  —  fi  <  (l/rmax)  [cf.  the  remark 
before  (42)]. 

Second,  we  derive  the  CRBs  for  real-valued  fi.  Let 


Vsi-  =  1/7  t]. 


(80) 


Define 


/  A 
F 1  = 


f4 


dp, 

dp 

dp, 

dr 


=  b 


=  j27r/i 


fiej2^hr 

f2pJ2^hT 


(81) 

(82) 
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